const char* dgemm_desc = "Naive, three-loop dgemm.";

/*#include <memory.h> */
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
/* This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values. */    
void square_dgemm (int n, double* restrict A, double* restrict B, double* restrict C)
{


    /* Reorganize A to be row-major */
    double* restrict Atemp;
    posix_memalign( (void**)&Atemp, 16, sizeof(double)*n*n);

    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            Atemp[i*n+j] = A[i+n*j];
        }
    }

    /* let's also copy stuff for B, so that it's aligned properly */
    //double* restrict Btemp = (double *)malloc(sizeof(double)*n*n); 
    
    double* restrict Btemp;
    posix_memalign( (void **)&Btemp, 16, sizeof(double)*n*n); 

    memcpy(Btemp, B, sizeof(double)*n*n);


    /* For each column j of Btemp */
    for (int j = 0; j < n; ++j) {

        int i;

        /* Vectorized part */

        /* For each row i, i+1 of Atemp */
        for (i = 0; i < n-2; i+=2)
        {
            /* Compute C(i,j), C(i+1,j) */
            double cij = C[i+j*n];
            double ci1j = C[i+1+j*n];
            double* restrict A_rowi = Atemp+(i*n);
            double* restrict A_rowi1 = A_rowi+n;
            double* restrict B_columnj = Btemp+(j*n);

            int k;
            for (k=0; k<n; k++) {
                cij += A_rowi[k] * B_columnj[k];
                ci1j += A_rowi1[k] * B_columnj[k];
            }

            C[i+j*n] = cij;
            C[i+1+j*n] = ci1j;
        }



        /* Non vectorized part */

        /* For each row i of Atemp */
        for (; i < n; ++i)
        {
            /* Compute C(i,j) */
            double cij = C[i+j*n];
            double* restrict A_rowi = Atemp+(i*n);
            double* restrict B_columnj = Btemp+(j*n);

            int k;
            for (k=0; k<n; k++) {
                cij += A_rowi[k] * B_columnj[k];
            }

            C[i+j*n] = cij;
        }


    }
    free(Atemp);
    free(Btemp);
}
